{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GPyOpt: Bayesian Optimization with fixed constraints\n",
    "\n",
    "### Written by Javier Gonzalez, University of Sheffield.\n",
    "\n",
    "## Reference Manual index\n",
    "\n",
    "*Last updated Friday, 11 March 2016.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we will learn how to solve optimization problems with fixed constraints. We will focus on problems where the goal is to find \n",
    "$$ x_{M} = \\arg \\min_{x \\in {\\mathcal X}} f(x) \\,\\, \\mbox{subject to}, $$\n",
    "\n",
    "$$c_1(x)\\leq 0 $$\n",
    "$$ \\dots $$\n",
    "$$c_m(x)\\leq 0 $$\n",
    "\n",
    "where $f: {\\mathcal X} \\to R$ be a L-Lipschitz  continuous function defined on a compact subset ${\\mathcal X} \\subseteq R^d$ and $c_1,\\dots,c_m$ are a series of known constraints that determine the feasible region of the problem. We will see the syntax that we need to use to solve this problems with Bayesian Optimization using GPyOpt. First we start loading GPyOpt and GPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import GPyOpt\n",
    "import GPy\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we will optimize the 2D Six-Hump Camel function (available in GPyOpt). We will assume that exact evaluations of the function are observed. The explicit form of the function is:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$f(x_1,x_2) =4x_1^2 – 2.1x_1^4 + x_1^6/3 + x_1x_2 – 4x_2^2 + 4x_2^4$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "func = GPyOpt.objective_examples.experiments2d.sixhumpcamel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imagine that we were optimizing the function in the intervals $(-1,1)\\times (-1.5,1.5)$. As usual, we can defined this box constraints as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "space =[{'name': 'var_1', 'type': 'continuous', 'domain': (-1,1)},\n",
    "        {'name': 'var_2', 'type': 'continuous', 'domain': (-1.5,1.5)}]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will be an standard case of optimizing the function in an hypercube. However in this case we are going to study how to solve optimization problems with arbitrary constraints. In particular, we consider the problem of finding the minimum of the function in the region defined by\n",
    "\n",
    "$$-x_2 - .5 + |x_1| -\\sqrt{1-x_1^2} \\leq 0 $$\n",
    "$$ x_2 + .5 + |x_1| -\\sqrt{1-x_1^2} \\leq 0 $$\n",
    "\n",
    "We can define these constraints as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "constraints = [{'name': 'constr_1', 'constraint': '-x[:,1] -.5 + abs(x[:,0]) - np.sqrt(1-x[:,0]**2)'},\n",
    "              {'name': 'constr_2', 'constraint': 'x[:,1] +.5 - abs(x[:,0]) - np.sqrt(1-x[:,0]**2)'}]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And create the feasible region od the problem by writting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "feasible_region = GPyOpt.Design_space(space = space, constraints = constraints)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's have a look to what we have. Let's make a plot of the feasible region and the function with the original box-constraints. Note that the function .indicator_constrains(X) takes value 1 if we are in the feasible region and 0 otherwise. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Grid of points to make the plots\n",
    "grid = 400\n",
    "bounds = feasible_region.get_continuous_bounds()\n",
    "X1 = np.linspace(bounds[0][0], bounds[0][1], grid)\n",
    "X2 = np.linspace(bounds[1][0], bounds[1][1], grid)\n",
    "x1, x2 = np.meshgrid(X1, X2)\n",
    "X = np.hstack((x1.reshape(grid*grid,1),x2.reshape(grid*grid,1)))\n",
    "\n",
    "## Check the points in the feasible region.\n",
    "masked_ind = feasible_region.indicator_constraints(X).reshape(grid,grid)\n",
    "masked_ind = np.ma.masked_where(masked_ind > 0.5, masked_ind)\n",
    "masked_ind[1,1]=1\n",
    "\n",
    "## Make the plots\n",
    "plt.figure(figsize=(14,6))\n",
    "\n",
    "# Feasible region\n",
    "plt.subplot(121)\n",
    "plt.contourf(X1, X2, masked_ind ,100, cmap= plt.cm.bone, alpha=1,origin ='lower')\n",
    "plt.text(-0.25,0,'FEASIBLE',size=20)\n",
    "plt.text(-0.3,1.1,'INFEASIBLE',size=20,color='white')\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.plot()\n",
    "plt.contourf(X1, X2, func.f(X).reshape(grid,grid),100, alpha=1,origin ='lower')\n",
    "plt.plot(np.array(func.min)[:,0], np.array(func.min)[:,1], 'r.', markersize=20, label=u'Minimum')\n",
    "plt.legend()\n",
    "plt.title('Six-Hump Camel function',size=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Six-Hump Camel function has two global minima. However, with the constraints that we are using, only one of the two is a valid one. We can see this by overlapping the two previous plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6.5,6))\n",
    "OB = plt.contourf(X1, X2, func.f(X).reshape(grid,grid),100,alpha=1)\n",
    "IN = plt.contourf(X1, X2, masked_ind ,100, cmap= plt.cm.bone, alpha=.5,origin ='lower')\n",
    "plt.text(-0.25,0,'FEASIBLE',size=20,color='white')\n",
    "plt.text(-0.3,1.1,'INFEASIBLE',size=20,color='white')\n",
    "plt.plot(np.array(func.min)[:,0], np.array(func.min)[:,1], 'r.', markersize=20, label=u'Minimum')\n",
    "plt.title('Six-Hump Camel with restrictions',size=20)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the modular iterface to solve this problem. We start by generating an random inital design of 5 points to start the optimization. We just need to do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# --- CHOOSE the intial design\n",
    "from numpy.random import seed # fixed seed\n",
    "seed(123456)\n",
    "\n",
    "initial_design = GPyOpt.experiment_design.initial_design('random', feasible_region, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Importantly, the points are always generated within the feasible region as we can check here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6.5,6))\n",
    "OB = plt.contourf(X1, X2, func.f(X).reshape(grid,grid),100,alpha=1)\n",
    "IN = plt.contourf(X1, X2, masked_ind ,100, cmap= plt.cm.bone, alpha=.5,origin ='lower')\n",
    "plt.text(-0.25,0,'FEASIBLE',size=20,color='white')\n",
    "plt.text(-0.3,1.1,'INFEASIBLE',size=20,color='white')\n",
    "plt.plot(np.array(func.min)[:,0], np.array(func.min)[:,1], 'r.', markersize=20, label=u'Minimum')\n",
    "plt.title('Six-Hump Camel with restrictions',size=20)\n",
    "plt.plot(initial_design[:,0],initial_design[:,1],'yx',label = 'Design')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we choose the rest of the objects that we need to run the optimization. We will use a Gaussian Process with parameters fitted using MLE and the Expected improvement. We use the default BFGS optimizer of the acquisition. Evaluations of the function are done sequentially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# --- CHOOSE the objective\n",
    "objective = GPyOpt.core.task.SingleObjective(func.f)\n",
    "\n",
    "# --- CHOOSE the model type\n",
    "model = GPyOpt.models.GPModel(exact_feval=True,optimize_restarts=10,verbose=False)\n",
    "\n",
    "# --- CHOOSE the acquisition optimizer\n",
    "aquisition_optimizer = GPyOpt.optimization.AcquisitionOptimizer(feasible_region)\n",
    "\n",
    "# --- CHOOSE the type of acquisition\n",
    "acquisition = GPyOpt.acquisitions.AcquisitionEI(model, feasible_region, optimizer=aquisition_optimizer)\n",
    "\n",
    "# --- CHOOSE a collection method\n",
    "evaluator = GPyOpt.core.evaluators.Sequential(acquisition)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we create the BO object to run the optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# BO object\n",
    "bo = GPyOpt.methods.ModularBayesianOptimization(model, feasible_region, objective, acquisition, evaluator, initial_design)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first run the optimization for 5 steps and check how the results looks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# --- Stop conditions\n",
    "max_time  = None \n",
    "max_iter  = 5\n",
    "tolerance = 1e-8     # distance between two consecutive observations  \n",
    "\n",
    "# Run the optimization                                                  \n",
    "bo.run_optimization(max_iter = max_iter, max_time = max_time, eps = tolerance, verbosity=False) \n",
    "bo.plot_acquisition()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See how the optimization is only done within the feasible region, out of it the value of the acquisition is zero, so no evaluation is selected in that region. We run 20 more iterations to see the acquisition and convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Run the optimization  \n",
    "max_iter  = 25\n",
    "bo.run_optimization(max_iter = max_iter, max_time = max_time, eps = tolerance, verbosity=False) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bo.plot_acquisition()\n",
    "bo.plot_convergence()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Best found value\n",
    "np.round(bo.x_opt,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# True min\n",
    "np.round(func.min[0],2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Done! problem solved within the fixed domain!"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
